Anharmonic effects in magnetoelastic chains 
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We describe a new mechanism leading to the formation of rational magnetization plateau phases, 
which is mainly due to the anharmonic spin-phonon coupling. This anharmonicity produces plateaux 
in the magnetization curve at unexpected values of the magnetization without explicit magnetic 
frustration in the Hamiltonian and without an explicit breaking of the translational symmetry. 
These plateau phases are accompanied by magneto-elastic deformations which are not present in 
the harmonic case. 

PACS numbers: 75.10.Jm, 75.10.Pq, 75.60.Ej 
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INTRODUCTION 

Coupling of electronic and elastic modes has been 
shown to play a crucial role in many condensed matter 
systems, most notably in the BCS theory of superconduc- 
tivity where the presence of the lattice degrees of freedom 
is crucial to explain pair formation [l[ . Another paradig- 
matic case is the so-called Peierls effect, where modula- 
tions in the charge or spin densities may appear due to 
the electron-phonon interactions (See e.g. |2| and refer- 
ences therein). More recently, phonon effects have been 
observed in many other strongly correlated systems, in 
particular in some magnetic systems which show plateaux 
in their magnetization curves 0] . 

Usually one expects to have an accurate description of 
an electron-phonon system by approximating the phonon 
potential with a quadratic function of the interatomic dis- 
tances between nearest neighbor ions on sites i and j, Sij. 
Within the same degree of accuracy, the dependence in 
of the hopping amplitudes and/or the magnetic ex- 
change constants is approximated as a linear function. 
This description works well in most of the cases since in- 
teratomic displacements are usually rather small as has 
been verified experimentally in many systems, like in the 
BCS superconductors. More recently however, a less con- 
ventional BCS superconductor, MgB2, has shown an un- 
usually high critical temperature, around 40°K, which 
could be the consequence of strong anharmonicitics both 
in the phonon potential and in the electron-phonon cou- 
pling 043- 

The relevance of anharmonic couplings has also been 
discussed in relation to a great variety of compounds, 
both from an experimental [8UlC| and a theoretical point 
of view (ill ] , including the family of pyrochlore oxide su- 
perconductors, AOS2O6 for A=Cs, Rb, and K ||, the 
heavy fermion superconductors PrOs4Sbi2, SmOs4Sbi2 



[9(, and some potentially thermoelectric materials such 
as XgGai6Ge3o (X=Eu, Sr, Ba) [l(|, etc. Another pos- 
sible relevance of anharmonicities is in the study of spin 
systems in high pulsed magnetic fields and Raman ex- 
periments (l2j . 

Apart from possible experimental motivations, the role 
of anharmonicities in the physics of low dimensional sys- 
tems is interesting in its own right and we investigate 
this issue in the present paper in one of the simplest and 
most paradigmatic one-dimensional systems: the XXZ 
Heisenberg chain. 

More precisely, we analyze in the present paper the 
effects of anharmonic (adiabatic) phonons in the spin- 
Peierls mechanism as well as the consequences on the 
magnetic properties of the XXZ Heisenberg chain cou- 
pled non-linearly to lattice deformations. The most im- 
portant consequence of the anharmonicity is that it pro- 
duces plateaux in the magnetization curve at unexpected 
values of the magnetization. For example a plateau 
at M = 1/3 of saturation magnetization appears with- 
out explicit magnetic frustration in the Hamiltonian and 
without an explicit breaking of the translational sym- 
metry [l3j], [lj]. Besides, magnetoelastic deformations 
appear in some particular cases with frequencies which 
halve that of the first harmonic, 2kp, as e.g. at M = 1/5 
(sec below). Similar conclusions should apply to more 
complicated models, since the effects of other interac- 
tions such as e.g. a next-nearest neighbor interaction 
would be simply to enlarge the extension of the plateaux 
phases and to modify the magnitude of the spin gaps [H[ . 



THE MODEL 

We start from the following spin-phonon Hamiltonian 
in the limit of large ionic mass M — > 00, the so-called 
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adiabatic limit 

H = J 5^(1 + A x fc + A 2 5 2 ) S t ■ S i+ i 

i 

-^E s >+E^) « 

Here <5i denotes the interatomic distance between site i 
and i + ft. is external magnetic field and S!j are spin 1/2 
operators. 

The dependence of the spin-phonon coupling on the 
interatomic distance Si has been expanded up to second 
order with coefficient A 2 . A Zeeman term is included to 
take into account magnetic field effects. 

The phonon potential energy in ([!]) is given by 



V(6i) = u 



(2) 



where 013 and cn take into account the anharmonicity of 
the interatomic potential energy. 

Generally, the properties due to the anharmonic os- 
cillations arise both from the addition of quartic terms 
in the potential energy and next-to leading terms in the 
spin-phonon coupling. In this letter we focus on the con- 
tribution of the anharmonicity in the spin-phonon cou- 
pling measured by A 2 ignoring the contribution of higher- 
order terms in the potential energy. We show below that 
it is the term quadratic in the lattice deformations in 
the interaction Hamiltonian that changes drastically the 
physics of the magnetoclastic XX Z chain. We expect 
that higher order terms in the potential energy (cubic 
and quartic) are inessential. 



BOSONIZATION DESCRIPTION 

Following the usual procedure in the low energy limit, 
we bosonize the spin degrees of freedom at fixed magne- 
tization M and the interaction term becomes 11511 



H sp _p h = / dx (A 1 5 M (x) + A 2 S M (x) 2 ) p(x) (3) 



where we have introduced the subscript M in 5m{x) to 
stress its dependence on the magnetization. Here p(x) is 
the continuum expression of the energy density 



p{x) = ad x <p + (3 cos(2kpx + V2ir(j>) + 



(4) 



where kp = ?(1 — M), a and j3 are constants and the 
ellipses indicate higher harmonics [16 1. 

The main contribution in the low energy limit comes 
from the constructive interference between the modula- 
tion term A\5m{x) + A 2 Sm(x) 2 and the most relevant 
part of p{x), i.e. cos(2kpx + y2~TT(f)). This operator 
has conformal dimension that depends on the Tomonaga- 
Luttinger parameter K(M, A)/2 where A measures the 



z-axis anisotropy in the XXZ model. Here we empha- 
size the dependence on the magnetization M and the 
anisotropy A. 

Let us propose a periodic pattern of deformations 
5m{x) with period L p , i.e. satisfying 5m (x + L p ) = 
Sm(x) (the lattice spacing a is set to 1 in what follows, 
so that L p is an integer). The most general Ansatz for 
the modulation term is given by 



Sm(x) = 53 $n(M) cos ( n- 

n=l 



e n {M) 



(5) 



where 6 n (M) are the amplitudes and 9 n (M) the phases 
of the different terms in the expansion. The upper sum 
index N w equals L p /2 if L p is even and (L p — l)/2 if 
it is odd. (In what follows the dependence of S n (M) 
and 9 n (M) on M is suppressed to ease the notation, i.e. 
6 n {M) -> S n and 6„ {M) -> 6 n ). 

From Eqs. Q and ([5]), we see that the product between 
the two terms is commensurate whenever the following 
relation is satisfied 



v 2n 
KF oc — , 



(0) 



which implies that the wavelengths of the modulations 
that could pin the relevant cosine term are related to the 
magnetization as 



Am 
1-M' 



(7) 



where 1 and m an arbitrary integer, the smallest 

possible that makes L p an integer. 

The Ansatz in Eq. (|5|) is verified a posteriori from the 
DMRG analysis, where it is seen that the modulation 
amplitudes S n and the phases 9 n depend strongly on the 
value of the magnetization M, some of them being zero 
in certain cases. 

Using this form for Sm(x), the interaction term J3]) 
takes the form 



H.«n- 



sp—ph 



2(AT„+1) 

E 

p=0 



X p / dx cos(p kpx + v 27r (j) + r p ) (8) 



where T p is a function of the phases 9 n in the expan- 
sion ([5]), and X p is a function of 6 n , 9 n and the coupling 
constants Ai and A 2 . 

This form of the interaction allows us to conclude that 
the spin Peierls effect takes place in the usual manner 
(see [l5J and references therein), since we have both the 
always commensurate term (p = in the above equation) 
cos(-\/27r <f>) and the 4fcp term, that provide together a 
dimerization of the lattice and a plateau at M = in the 
magnetization curve. 
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M = 1/5 


M = 1/3 


M = 1/2 




5 


6 


8 


N w 


2 


3 


4 


z±2 


5,10,... 


6, 12,... 


8,16,... 


z 


3 


4 


6 


possible 
frequencies 


kF and 2kF 


2k F 
or 

kF and 3k p 


2kF and 4A;f 



TABLE I. Possible frequencies for the lattice deformations for 
magnetizations M = 1/5, 1/3 and 1/2, obtained from the Eq. 



For finite magnetization, using Eqs. ©-© and us- 
ing the commensurability condition that arises from 
p kF /2tt E Z, one obtains the following condition for the 
frequencies in (JS|) to pin a relevant perturbation 



(z ± 2) (1 — M) = 4 x integer 



(9) 



where z is an integer that runs through all the frequencies 
that appear in the lattice deformation Eq. and its 
square, i.e. z = 0, ...,2N W . In Table [J we show some 
examples that we analyze in what follows using DMRG. 

One should stress that in the present case the situa- 
tion is rather different than in previous studies of spin 
systems in a magnetic field, such as in the case of spin 
ladders, magnetoelastic zig-zag chains, etc., since now 
the perturbing operator that would be responsible for the 
plateau is relevant, independently of the values of the mi- 
croscopic parameters. This may seem to imply that con- 
dition (j9]) is also sufficient, but bosonization alone does 
not provide the actual values of the amplitudes of the 
different Fourier components of the deformation we pro- 
posed and it remains to be checked that they are indeed 
non vanishing. In order to answer this question we need 
to use DMRG as we describe below. 

From the above analysis, we predict that the magne- 
tization curve may present new features related to the 
frequencies which appear in the Fourier decomposition 
of the elastic deformation <j5j> for some given values of 
M, such as M = 1/5, M = 1/3, M = 1/2, etc. Since 
these frequencies pin the very relevant term cos(v / 27r (f>), 
plateaux at these values of M are expected to show up 
even for a small anharmonicity A2- In such cases, the 
plateaux widths Gap(M 7 A) should scale as 



Gap(M, A) oc a 1 /(2-^,a)) 



(10) 



where A is the coupling constant associated to the rel- 
evant cosine term in H sp - p h and d(M, A) is the scal- 
ing dimension which can be computed from the Bethe 
Ansatz solution [17| . The coupling A is a function of the 
anharmonic amplitude Ai and its functional dependence 
though not predicted by bosonization, can be computed 



numerically as we show below. From now on we will 
concentrate in the isotropic case A = 1. 

DMRG ANALYSIS 

This is the general setting obtained from bosonization 
which provides the qualitative picture expected when an- 
harmonic effects play a role. To have a complete and 
more quantitative picture we study the system using ex- 
tensive DMRG computations. More specifically, we com- 
pute the ground state energy E(S^ otal , h = 0) of Eq. ([T]) in 
the complete set of S^ otal subspaces using periodic bound- 
ary conditions, and keeping just 300 states it was shown 
to be enough to assure the accuracy of the calculation. 
As usual, adding the Zeeman term, we solve the equa- 
tion E(Sf otal , h) — E(Sf otal +l,h) to obtain the normal- 
ized magnetization M ~ 2S Z /N where the plateaux are 
showing up. This procedure allows us to compute the 
actual width of the plateaux and their scaling behavior, 
the deformation patterns and fractional excitations for 
the different plateaux. 

Let us analyze in detail the situation at M = 1/3, 
where we expect to have a plateau. In this case kp = 7r/3 
and our Ansatz for the modulation ([5]) takes the form 

S 1 / 3 (x) = Si cos(fcp x + 6{) + 82 cos(2kp x + 82) 

+ 83 cos(3fci? x + 3 ) (11) 

which leads to the perturbation Hamiltonian 



H 



sp—ph 



A1/3 cos (\ 



2tt( 



(12) 



where X1/3 and T 1 / 3 depend on Ao and Xq which are the 
only two commensurate terms in Eq. ([S]) at magnetization 
AI = 1/3 (see Fig. []}. The dots in the equation above 
indicate less relevant terms, which can be safely discarded 
in the presence of the more relevant term oc cos [V^tt <£>) . 

The couplings appearing in (|12[) have a lengthy ex- 
pression in terms of the strengths of the spin-phonon 
couplings Ai and A 2 but also on the <5„'s and on the 
relative phases 6*„'s, whose values cannot be extracted 
from the bosonization analysis alone. To further proceed 
we now resort to the numerical analysis of the system us- 
ing DMRG on large systems which allows us to estimate 
all these parameters in a self-consistent way. 

The lattice deformations can be calculated in a self 
consistent way. Minimizing the ground state energy and 
imposing the following constraint 







we obtain 
JA 1 



T, k (Sk-S k+1 ) (u> +2JA 2 (S k -S k+1 ))- 
E fc ("o+2JA 2 <S fc -S fc + 1 >)- 1 



(13) 



— (Si ■ Si+i) 



(uj + 2JA 2 (S 1 - S t+1 )) 



(14) 



4 




FIG. 1. M vs. h for Ai = 0.6, A 2 = 0.4 and different system 
sizes. (N = 30,36,48 and 60). The bold purple line cor- 
responds to the extrapolation to the thermodynamic limit. 
The plateaux at M = and 1/3 are clearly observed, while 
for M = 1/5 and 1/2, it is hard to conclude if they survive in 
the thermodynamic limit. Note that for N = 30, M = 1/5, 
1/2 are not commensurate. The inset shows the finite size 
scaling of the width of the plateau at M = 1/3. Its finite size 
scaling is expected to follow width(N) = width(oo) + AN~ B . 



We start from an arbitrary chosen initial set of deforma- 
tions {<5^ ^} to be varied and determined self consistently. 
For a given set {o^- 1 } we determine the corresponding 
ground state and then we compute a new set {<5^ Ar+1 ^} us- 
ing (|14j) which we use again in the Hamiltonian. Iterating 
this procedure, we finally obtain a fixed point configura- 
tion of the deformations S^ N+1) ({S^}) = sf\ 

From the DMRG data, we observe that for M = 1/3 
only the 2fcp mode contributes to the lattice deformations 
(see Fig. [2]), so that we can safely set S\ = S3 = 0. As 
for the phase 2 , it is negligible within the numerical 
precision so we set it to zero in what follows. With this 
input from DMRG, we get the following expressions for 
the bosonization parameters, i.e. for the amplitude \\/3 
and phase r 1 / 3 in Eq. (|T2|) . 



A1/3 oc 



Ax 8 2 



A 2 S\ 



Tl/3 = -7T/3 



(15) 



Here a word is in order: To analyze the scaling of the 
gap we need to identify the effective coupling constant 
associated to the perturbing operator responsible for the 
opening of the gap. Since the term proportional to A\ is 
present for all magnetizations, it does not play a role in 
the gap opening and we can then identify the coupling 
constant governing the scaling of the gap in (fTOj) as A a 
A 2 Si 

On the other hand, we can extract the deformation 
amplitude as a function of A 2 from the numerical data, 
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= 0.3 




A 2 


= 0.4 



1 2 

w n /k F 



FIG. 2. Amplitudes 5 n (M) (see Eq. ©)) as a function of the 
frequency w n — 2nn/L v in units of &f (with Ai — 0.6) for 
M — 1/2, 1/3 and 1/5. The peaks indicate which frequencies 
contribute to the deformation pattern. The 2&f peak is al- 
ways bigger because to linear order in 5m(x) it contributes to 
the energy for all magnetizations. 



which after a finite size scaling analysis and a square fit 
leads to S 2 = a + bA 2 + cA% with a = 0.110, b = 0.098 
and c = 0.551. Now that we have the dependence of 
the effective coupling A1/3 on the anharmonicity A 2 we 
can analyze the scaling of the spin gap (the width of the 
plateau), which should scale as in (fT0|) . 



log Gap = a logAi/3 + cte 



-2.5 r 
-3.0 
log Gap— 3.5 
-4.0 
-4.5 L 




Bosonization 
DMRG 

a = 0.778 

reference 
curve 

a = 0.77 



FIG. 3. Logarithmic plot of the Gap(\): The blue dots cor- 
respond to the reference curve Gap(X) = A 77 with the gap 
obtained from DMRG, while the red dots correspond to the 
gap obtained from DMRG vs. the values of A extracted from 
bosonization. The value 0.77 is obtained from the Bethe 
Ansatz solution 

In order to compare both approaches, we need to use 
the relation (jT5"j) between A1/3 and A 2 , together with the 
values of S 2 obtained form DMRG. Following this ap- 
proach, in Fig. [3] we show a logarithmic plot of the gap 



5 



vs A1/3 using the values of A1/3 obtained by bosonization 
and those of the gap by DMRG (red points). We show a 
linear fit to obtain the exponent in eq. (|10p and compare 
with a reference line (blue points) to show the agreement 
of both approaches. 



CONCLUSIONS 

In conclusion, we have described a new mecha- 
nism leading to the formation of rational magnetization 
plateau phases, which is mainly due to the anharmonic 
spin-phonon coupling. We have shown that its role is to 
pin magneto-elastic deformations which are not present 
in the harmonic case. By means of bosonization we have 
shown that the inclusion of the anharmonic spin-phonon 
coupling gives as a contribution a relevant operator that 
is responsible for the plateau in the magnetization curve 
for certain commensurate values of the magnetization 
M. We have performed extensive DMRG computations 
to complement the analytical computations, since the 
bosonization approach alone does not provide the ac- 
tual values of the amplitudes of the different Fourier 
components of the lattice deformations. In particular, 
we have analyzed in detail the situation at M = 1/3 
where we have computed the plateau width as a func- 
tion of the anharmonic coupling, to extract the scaling 
dimension of the relevant operator that opens the gap. 
Finally, we have seen that the exponent obtained from 
the DMRG computations and the one obtained from 
the Bcthe Ansatz through bosonization, are in excellent 
agreement, providing further support to our results. 
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